
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#COMPUTE ELASTICITY η(SP,BE) FOR DIFFERENT ALPHAS
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------

#Load functions for simulations
#------------------------------------------------------------------------------
include("../01_Baseline_Model/01_Model_Structure.jl");
include("../01_Baseline_Model/02_Functions.jl");
include("../01_Baseline_Model/10_Simulations_Main_fx.jl");
include("../01_Baseline_Model/11_Simulations_Moments_Unconditional.jl")
include("../01_Baseline_Model/12_Simulations_Wrapper.jl")

using Interpolations, Optim, Roots, JLD, Plots
using Distributions,  StatsBase
#------------------------------------------------------------------------------


#Set options
#------------------------------------------------------------------------------
name_file        = string("../01_Baseline_Model/model_data/CE_model_structure_alpha_",0.028,".jld");
ce               = load(string(name_file), "ce");
avec             = [-0.034; -0.032; -0.030; -0.029; -0.028; -0.027; -0.026; -0.024; -0.022];
elastAVG_alpha   = zeros(length(avec));
T                = 200_000
N                = 1
#------------------------------------------------------------------------------


#------------------------------------------------------------------------------
#Loop across α
#------------------------------------------------------------------------------
for jj=1:length(avec)
    println(avec[jj])
    #Load Economy
    name_file = string("../01_Baseline_Model/model_data/CE_model_structure_alpha_",-avec[jj],".jld");
    ce = load(string(name_file), "ce");

    #Simulate the economy once to compute some moments
    N = 1
    y_sim, type_sim, b_sim, π_sim, def_sim, def_pol, def_id, q_sim, SP_sim, SP_HR_sim, c_sim, tb_sim, m_sim,m_sim_ctfl, ζ_sim,
    dlnSP_sim, dBE_sim, b_sim_tmp, inf_bp =
    simul_econ(ce,200_000,N, T0=5_000);
    RP_sim = SP_sim - SP_HR_sim;
    π_std      = std(π_sim[(type_sim.==2) .* (def_sim.==0)])
    ζ_avg_C    = mean(ζ_sim[(type_sim.==1).* (def_sim.==0)])
    ζ_std_C    = std(ζ_sim[(type_sim.==1) .* (def_sim.==0)])
    ζC_l       = ζ_avg_C-ζ_std_C
    ζC_h       = ζ_avg_C+ζ_std_C
    #Simulate the Economy and compute the elasticity
    y_sim, type_sim, b_sim, π_sim, def_sim, def_pol, def_id, q_sim, SP_sim, SP_HR_sim, c_sim, tb_sim, m_sim, m_sim_ctfl, ζ_sim, dlnSP_sim, dBE_sim, b_sim_tmp, inf_bp,
    elasticity, ΔSP_sim,ΔBE_sim,Δπ_pol = 
    simul_econ(ce,T,N, T0=5000, π_shock=π_std)
    #Filter observations and compute the average elasticity
    filter_sp           = filter_fx(def_sim,type_sim,ζ_sim, SP_sim, ζ_min=ζC_l, ζ_max=ζC_h)
    elastAVG_alpha[jj]  = mean(elasticity[filter_sp.==1])
end
#------------------------------------------------------------------------------


#Create Plot
#------------------------------------------------------------------------------
Plots.plot(avec, elastAVG_alpha,  xlabel="\\alpha", ylabel="Elasticity", linestyle=:dash, linewidth = 2, label=(""), color=[:black],  markershape=:circle, markersize=2.5, markercolor=[:black], markerstrokealpha=0.4)
Plots.savefig("figures/elasticity_alpha_comparisons.pdf")
